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(57) Abstract: Movement of material in an organism, such 
as a drug injected into a brain, is modeled by a uniformly 
structured field of static constants governing transport by 
moving fluid and diffusion within the fluid. This supports 
planning of material introduction, (e.g., infusion, perfusion, 
remjperfusion, injections, etc.) to achieve a desired distri- 
bution of the material, continuing real-time feedback as to 
whether imaged material is moving as planned and will be 
distributed as desired, and real-time plan modification to 
improve results. 
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A METHOD AMD APPARATUS FOR TARGETING MATERIAL DELIVERY TO 
TISSUE 



BACKGROUND OF THE INVENTIOIV 
The present invention relates to the delivery of material to tissue within a subject or 
patient. In particular, deUveiy of materials such as therapeutic, image enhancing, bio-active, 
pharmacological, nanotechnical or otherwise active materials is enhanced, particularly under 
real-time observable imaging systems (such as MRI, sonograms, and X-ray fluoroscopy). 
The invention forther relates to the field of predictive mass transport or diffiision analysis for 
enhancing material delivery within a subject or patient. 

It is increasingly common to administer a drug or other material to a care&Uy targeted 
part of the body, rather than to insert it into the bloodstream and rely on some of it finding the 
target by carriage through the circulatory system. This targeted delivery has multiple 
advantages, when it can be performed effectively. Among the advantages of targeted 
deKveiy are that much less of the drug is needed, which itself represents a double gain. The 
drug itself is often costly, so reducing the volume of drug used in a treatment can represent a 
significant cost savings. It is also rare that any drug is wholly without negative effects ('side 
effects' with reference to the desired result), mere these adverse effects arise only where 
the drug reaches a specific non-targeted tissue, restincting tiie drug to a target that does not 
include that non-targeted tissue avoids the side effects completely. Even where this complete 
exclusion is impossible, the side effects may be far more acceptable if limited to a smaU 
region around the target, with reduced impact on the body at large, while delivering the 
desired result on the target tissue at full strength. Even where the side effects are not directly 
life-threatening, it can be important to avoid them: for example, to keep cancer chemotherapy 
drugs fi-om the sites where they cause nausea and hair loss is good both for patient morale 
and for patient persistence in taking the drug. The importance of targeted delivery wiU 
increase for nanodevices (that is, devices with dimensions that are measured with one or two 
orders of magnitude of nanometers), which will often be designed for highly specific activity 
in a particular environment. Apart fi-om the waste of resources in feilure to reach the target 
site, their action in unintended (and less studied) sites may be hard to predict. 
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However, this format of targeted material delivery adds to the traditional complexity 
of computing dosage and deUvery rates. Instead.of a single figure of blood concentration. 
controUed by the rates at which the drug enters the circulatory system and at which it 
diffiises, flows and is absorbed, metabolized or excreted by various tissues, concentration 
5 becomes a distinct time- 

varying value (number) at each part of the body.. Since the effects of a drug vary in complex 
ways depending upon local levels of concentration (scopolamine, for example, is an anti- 
nausea drug in a narrow range of levels), its administration must ensure that at active points 
10 in the target tissue the concentration is correct acrossthe entire region of targeted tissue to 
obtain the desired effect, while at other non-active treatment points minimizing undesired 
effects (often, but not always, by minimizing concentration at such points). Planning deUvery 
commensurate with desired treatment effects thus requires a difficult prediction task. 
The need for such prediction of delivery and diffusion profiles has increased, 
15 particularly with the recent development of methods of tracking a diffiising or flowing 

material in real time in the patient (see U.S. Patent No. 6.026.316) and with the advent of 
direct drug infosion techniques (see U.S. Patent No. 5.720.720. which describes a catheter- 
based technique for high^flow microinfosion, and U.S. Patent No. 5,735,814, which discusses 
drug infusion into braintissueby means of an implantable pump and catheter). The 
20 physician can thus observe the changes in: concentration in the various tissues around the 
entry point and delivery region, and modify plans.a<3cording:to.observed events. It is 
important to note that the abiUty to modify or alter thesignificant results are not effectable 
instantaneously. Unlike an artist applying paint to a canvas, where the result is immediate, 
the physician must control the administration process according to events that wiU 
25 subsequently provide observable or therapeutic effects that result over a period of seconds or 
minutes. However, the decision itself must be immediate, so the prediction of consequences 
must be available immediately in advance of the alteration of procedures. 

The physician's brain, or a computer in assistance, must model the concentration 
dynamics prospectivdy much faster than they occur in real time, to be usefiil in real-time 
30 decision making. However, current computational methods of predicting concentration 

dyaamics in tissue take longer than the actual events in the process of delivery. These current 
methods rely on two steps, both of which are slow. 
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The first step, once a scan of the target region is available (in either preoperative or 
real time during initial steps of the.medical procedure), is to build a structural model of the 
tissue structures present. This requires first the labeling of points according to the type of 
tissue present (this process is often caUed 'segmentation', since it categorizes the points into 
5 three-dimensional 'segments'). Figure 1 iUustrates this with a 2D sKce of a brain scan, with 
the Globus Pallidus Medialis on each side extracted and marked visuaUy (101) by a uniform 
gray shade. Figure 2 shows the layered 3D region constructed fi-om such slices. The planner 
then creates a best-ft geometric model of each structure present, such as bone, hippocampus, 
cortex, etc. Figure 3 shows such a 3D model for one of the segments identified in Figure 1. 
10 Note that 'model' in this sense is a description of a geometrical shape, by (for instance) 

specifying vertices and feces, rather than a statistical model of a relationship, found by such 
methods as least squares. The decision as to what model fits the data best is somewhat 
heuristic: usuaUy the objectives in fitting a geometric surfece model to a bone are that points 
inside it should mostly be 'probably bone' on the evidence of local scan values, that they 
1 5 should form a connected region, that points immediately outside should be 'probably not 
bone', and that the surfece should be reasonably smooth. (This last criterion tends both to 
reduce the impact of noisy data, and to aUow a model that uses fewer vertices and feces.) 
Methods for constructing such a surfece model range from local definition of a surfece that 
separates points according to whether they are above or below a threshold value, such as the 
20 Marching Cubes technique [W E Lorensen and H E Cline, System and method for the display 
of surfece structiwes contained witliin the interior region of a solid body, US Patent No. 
4,710,876] to active 'baUoons' that move over the 3D image and attach themselves to 
boundary-like points, while resisting extremes of bending, (See for example L Cohen, L D 
Cohen, and N Ayache, "Using deformable surfeces to segment 3-D images and infer 
25 -differential structures," CVGIP: 19, Image Understanding 56(2):242-263, Sept. 1992 ) 

These models are divided into finite elements with simple geometric forms, such as 
tetrahedra (as in Figure 3) or skewed cuboids, spheres or other geometric or mathematical 
shapes. On each such element, a partial differential equation governing concentrarion 
dynamics, which by the definition of 'differential' involves values at an infinitude of points, 
30 is approximated. This approximation is by a system of equations with a smaU or at least 
controlled number of variables. Typically each variable multiplies a fixed fiinction of 
position before it is added to an approximation of the concentration fimction. In the simplest 
cases of current art, these fimctfons may be constants and linear fimctions, such as fimctions 
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proportional to X. proportional to and proportional to^ axes in the volume. In more 
complex approximations, they may be polynomial fimctions, wavelets, etc. The interactions 
between the scalar coefficients of these fimctions. within a finite element and between 
neighbors, are chosen to reflect the local rates of diflBision. flow, absorption, metabolism, 
excretion, etc., in a 'lumped' feshion over each finite element, so that the second step, 
computing the evolution of the approximated concentration fimction. gives results 
approximating those that would be predicted by a solution of the underlying partial 
differential equation. 

The model-construction step is a relatively long procedure. It is not easUy automated, 
and yet it must be complete before any predictions are possible. (Figure 3 represents many 
hours of technical work.) Currently it takes much more time than would be acceptable in 
clinical practice, creating a delay between a preliminary scan and the administration of a 
drug. The scarce and costly expert labor required for performing this step is another 
obstruction to deployment. When the prediction is complete, the predictwe calculations may 
stUl remain too slow for real-time guidance of the drug administration process. 



Brief Summarv of the I nvention 
The methods described in U.S. Patent No. 6,026,316; and instruments used mU.S. 
Patent No. 5,964,705 and copending U.S. Patent Applications bearing attorneys' docket 
20 numbers 00723.026US1, 723.029US1, and 723.030US1 filed on March 21, 2000, mal^e it 
possible to administer a trackable material at one xDr more chosen points in the organis^, and 
record the resulting change over time of the material's concentrations at a large numb^ of 
points. From this information may be estimated the characteristic rates at which the tissue at 
these points permits the diffusion and other processes by which the concentration evolves. 
25 For example, if the concentrations at two neighboring points are always very similar, 

diffiision rates between those points must be larger than diffiision rates between nearby points 
where substantial differences in concentration persist. In quantifying these deductions or 
estimates, one must treat all the terms of the concentration dynamics simultaneously, since 
the action of each term of the concentration has had impact on the available data. These 
30 estimates may be refined by certain fonns ofa priori knowledge; for example, bone can be 
recognized by static aspects of the scan data, without reference to recorded concentrations, 
and diffusion and flow are known to be slow in bone. Using such scan data, with 
segmentation infomiation where convenient, but without constructing a geometrical model 
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specffically representing the segmented tissues, we construct a field of coefficients (such as 
directioD-dqiendent diflSision rates, flow velocity of any fluid which transports the material, 
absorption rates, excretion rates, and any other processes important to the concentration 
dynamics). This creates a highly regular structure (e.g.. Figure 4, and Figure 5) of data and 
evohition rules, varying from point to point only by the static local coefficients. We then 
model the evolution of the concentration either by a field of values (where the static local 
coefficients control the disappearance rate of ttie material and the interactions between 
neigjiboring concentration values), or by a population of moving particles whose probability 
of transfer between neighboring points, or of disappearance, is governed by static local 
coefficients. Figure 6 illustrates this standard method for the case of a 2D field evolving by 
diflfiision only. A group of modeled particles startmg at point 601 will spread out more 
rapidly into a dilute cloud than a group starting at point 602, if their motion is computed 
according to the constants in Figure 5. 

The creation of such data structures may be largely automated, foUov/ing the flow 
indicated in Figure 7. The process is computation-intensive, but can be greatly speeded by the 
use of parallel computing: an option is to use a highly pipelined machine such as a Cray 
mainfi^e computer, but — particularly with the rapid increase in machine speeds — a 
multiprocessor PC will suffice for an increasing portion of the fonctionaJity of the present 
invention. For a small fi^ion of the cost of the imaging device such as an MRI scanner that 
is providing the data, a chip or chipset optimized for this particular task will enhance the 
power of the images in guiding intervention. The automated process reduces delay and costs 
to levels acceptable for clinical deployment, and the concentration computations are of a type 
highly susceptible to computational speedup by parallelization methods, such as the use of 
pipelined vector architecture for a field of concentration values, or of 'single instruction 
multiple data' (SIMD) architecture to track a population of particles. This process enables 
computation to predict the results of particular administration strategies fest enough so that 
plans can be compared. An optimal plan can be chosen, in a clinically acceptable time. The 
plan can be implemented while the material administration procedure is in progress to track 
deviations Grom the expected evolution of concentrations, to deduce corrections in the static 
coefficients, and to offer both revised predictions of the results of the current strategy and 
predictions of possible alternative strategies. All of this may be done in time to be usefiil to 
the controller. 



PCT/USOl/14633 

WO 01/852^J 



Figurel. A2Dscanimage,representingthe3DdataactuaUyavailable,mw 
particular tissue (the Globus Pallidus Medialis or GPm, shown uniform gray) has been 

identified. ,n 
5 Figure 2. The 3D form of the GPm in the same patient, buih 6om shoes like those m 

Figure 1, as an unstructured volume. 

Figure 3 The 3D region shown in Figure 1. dwided into geometric elements on 
whichconcentrationvariationmaybeapproximatedbyaweightedsumofasmallsetof 

standard functions. 

10 Figure4. A regular 2D grid of scan samples, lepresenting the 3D analogue, each of 

which is shown with a gray shade corresponding to its classification as to tissue type. 

Figure 5. Aregular 2D grid of points, representing the 3D analogue. To each point 
are assigned coefficients corresponding to its type (as shown in Figure 4) for the various 
processes by wWch the materialofinterestistransported or is removed from the ^st^^ 

15 simplicity, a single coefficient is shown. 

Figure 6. A randomly-moving particle at point 601, where according to Figure 5 the 
diffiasion constant is high, is modeled as more likely to move to a neighboring point than a 

particle at point 602, where the constant is tow. ' 

Figure?. The flow chart ofthe computation and interaction procedure of the 

20 invention. ' ^ ■ ■ - 

In the practice of the present mventfonj a^giot. of the body b treated as a fieU f^f 
local propertie. such as difiusion and abson-Uon constai«s, rather than as an assembly of 
organs Wc follow the convention that a 'field- is a specification, at each point in a domam, 
23 of a scalar, vector or tensor quantity characterizing a property, or of several such: the value of 
P at a point . is the ordered list oflocal property vataes at This field vary in ttae (ftr 
stance, with heart or lung motion) as weU as in space. Afficted by this field, but taken as 
without influence on it, is a velocity field v of moving fluid. TMs moving fluid may in some 
oases be a steady or periodic motion due to bodUy processes, but in many cases of interest the 
30 velocity field is more often a transient, time-varying flow driven by a pressure difference 
between insertion devices, concentration variations, osmotic pressure, and backgmund 
pressure in the organism. At some stage after insertioh ceases or release ofthe material 
ceases fte flow velocity v may become negligible. Fimdly. one object of interest .s the 
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concentration field C of the concentration. of a dififijsing material such as (but not restricted 
to) molecules, particles, globules, a dispersed or suspended phase, cells or microscopic 
devices, dissolved or in suspension in themoving or static fluid. In some cases of interest for 
the invention, this field varies in time. In one implementation, we treat the concentration as 
influenced by, but without influence upon, the fields P of local properties and v of flow 
velocity. One method of the invention is (a) to estimate P, (b) to estimate v, and (c) to use 
these estimates to predict the evolution of C, and (d) to use this capability of prediction for 
various purposes of medical or biological utility. In certain cases, an impact of C upon P or v 
(for instance, where a difiusing drug stimulates a change in the mechanical properties of 
tissue) requires an iterative process of improving the estimate (a) or (b) in the light of the 
results of (c), using the result for a repeat of (c), and so on. 

The estimation (a) ofP may be done in any manner. At least two wa5's will be 
specifically described. A first method is to use scan data to classify ('segment') each point in 
space as belongmg to the specific tissue tyije for which the scan value found there is 
characteristic, and to assign corresponding coefiBcient values based on pre-established 
knowledge of that type of tissue. Such assignment is iUustrated (in 2D analogue) in the 
transition from Figure 4, which shows points classified by tissue type and visualized by 
distinct gray levels (or color levels or vector symbols), to the data structure in Figure 5, which 
shows difiiision constants associated with the tissue type at each point. (In general, more 
than the single number shown is necessary to encode the various parameter forming part of 
P) ........... 

In some cases, the segmentation process can be improved by smoothing or otherwise 
processing the raw data or mathematical models, or even by constructing best-fit geometric 
models of the segments or such boundaries as the surfece where tissue meets a fluid-filled 
ventricle, to enforce agreement with such a priori knowledge about a segment's shape as that 
the outer membrane of the cerebral cortex is a folded surface without holes and without tubes 
connecting parts of it, but such geometric models are incidental to the present invention. With 
the current low speed and degree of automation of such models, they do not at present form 
part of a preferred implementation. 

A second example of a way to estimate i* is to inject material {e.g., scout or testing 
material, MRI fieU modifying, sonogram enhancing, inert materials, etc.) whose 
concentration is detectable in a scan. The fieW of densities is then detected as a fonction of 
position and time, and the properties computed directly from the changing values detected. 
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Since these changes in G are affected also by non^negligiDle values of v, we describe 
techniques for such computation in a later section. : - ■ 

Similarly, for step (b), v may be estimated in several ways. One way is, agam, to 
begin with a segmentation, and to provide an equation assumed or estimated to govern the 
5 flow withtermsintheequationforthee&ctsofthetissueitisflowingthrough. Withthe 
assumption of an incompressible fluid (usually appropriate when a diffiising drug is mjected, 
though less exact for lung air through which an anesthetic dififoses), fluid transport is 
governed by 
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where at a particular point x of the tissue the value ^x) is Ihe pore S^pOc r) « the 
(spatially a.d temporally varying) fluid pressure. is flie tensor of hydmuBc pem«abi>ty, 
and uix) is the spatial displacement of that point of the tissue under Ae forces fiom the flow. 

The quantities v. « and p are solved for«.llectively. usmg boundary conditions such 
as that fer from the injection site, p has Us bacl«round value for the organism. Equation (1) 
cannot of course be soh,ed alone: fer quasi-steady flotv (negligible acceleration) we have the 
force balance equation 

. . V.(A.e)-V(#) = 0 (2) 

where A is the tensor of constitutive moduU for the tissue and 



i 








2 









an 



is the tissue strain tensor, af the pore fraction ^ is modified by change in pressure, 
equation for this influence on equation (2) must be added.) We require also the effective 
continuity equation 



25 V.W + V- 



a, 

where Q is determined by the net rate F. of fluid loss from blood vessels into the tissue per 
unit tissue volume, and the absorptionrateF.offluid by lymphatics per mutvolumeoft.^^^^ 

Numerous modifications to these equations will be apparent to those skilled m the art. 
Mathematically modified, situationally modified, structuraUy modified, compositionally 
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modified, environmentaUy modified or otherwise fiinctionally modified equations or 
algorithms can be used in the practice of this invention, and these modified equations are 
assumed and intended to be included within the recitation and descr5)tion of these equations 
or algorithms except where the tenn "consisting of directly precedes the equation or 
algorithm. For example, if tlie.diflEusing substance is physiologicaUy active, it may modify 
the pore fiaction 4> (and hence the hydraulic penneability K) by causing cells to change shape, 
or the fluid loss rate Fy by constricting or relaxing the blood vessels. If the tissue is 
sufficiently rigid, we may neglect the force balance equation (2) and reduce transport to 

A potential result of this invention is not to give a universal system of equations 
governing transport in tissue, or indeed a universal approach to solving such equations, but to 
explicate a method that integrates such equations and solution methods into a system jfor 
interactive simulation. We will illustrate the method using equations (3) and (4). In this case 
the field of properties P would be the coUective name for {ip, Q, K). In other instances it 
would include (for instance) the tensor A of constitutive moduli for the tissue, the local 
coefficients governing the influence of the concentration Con pore fi-action and so forth. 

The values of fields such as 25: and <^ (or. if necessary, their dependency laws on other 
quantities in the equations) can typicaUy be taken fi-om pre-established knowledge of the 
properties of different tissue types, though these values may also be derived from 'scout' 
injections or revised in the hght of incoming data during a procedure. The equations may be 
solved to find V by a variety of (mostly approximate) numerical methods, given boundary 
conditions specified by a plan to inject material at a given forced rate of flow, perhaps time 
varying, or at a given pressure diflferential, perhaps time varying, over the bacl^ound 
mtemal pressure of the body. Equations may also be influenced by specific blood fectors and 
properties and conditions, and data from blood tests or other physiological tests may be used 
to modify the equations. The tests may be performed before or during the procedure, and 
certain tests may be performed outside of the patient for determining modification fimctions 
or modified constants based upon those blood properties. 

The equations may also be modified for convenience in approximate sohztion For example, 
the partial differentials dldK.^dldx^ and didx, applied to each quantity ^ (standing for «! v, 
etc., as appropriate) may be replaced by the i^specthAe partial differ^ce operator 
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(,(x.^.rep.x,,X3)-.U,x,,X3))/W, (.U 
{q{x„X2,x,+step)-q{xi,X2,xM^tep], 

and time derivatives sinularly replaced by time difference operators. This creates a system of 
equations involving only values at a finite set of grid points at a finite set of times. Such a 
; system5maybesolvedl^(forinstance)takinganinitialestimate,appro^^^^^ 

for the value of each adjustable quantity at each grid point, computing the nearby partial 
differences, assessing the Mure to satisfy S at that point, adjusting each scalar m the 
directionofsatisfying^better.andrepeatingtheprocessfortheadjustedvalues. (The 

amount by which to change the values at each step, for convergence to a solution and rapid 
0 convergencetoasolution,isthesubjectofalargesectionoftheN«mericalAnalysis 

Uterature.) Altematively. the region in which the system S is to be solved may be subdivided 
into smaller regions, on each of which each quantity q is approximated by a hnear 
combination c.g..c,,,...- + c...ofafinitecoilectionofMstandardZ,..^^^^^^^^ 

which may be polynomials, wavelets, or other mathematicaUy usefUl forms: integrating the 
system S over a sub-region gives a system ^ of equations relating the unknown coefficients 
(c.,c,....,c^), rather than the point values of the quantities q. (Solution for each subregion, 
on wWch ihe fimctions of interest in the differential equation require an infinity of values for 
fun specification, is thus reduced to determining a finite set of coefficients. Such techmques 
^t^crcforco^c^finiteelementmethods. The current drug diffiision modeling 
techniques, descrfted above, depend on detertmhing anatomy and dividing particular tissues 
into finite elements, like the tetrahedra shovm in Figure 3. It is an essential feature of the 
present invention that we do not depend on this, though we may use automated methods to fit 
a finite element decomposition to the continuous fields of quantities The systemgmay 
then be soh^ed by iterative adjustment fi-om an initial guess, or otherwise. 
15 A second approach to performing step (b) is to inject comparable fluid, which is 

visible in some way to a scan (as by MRI, sonogram. X-ray fluoroscopy, etc.), track its 
motion without assmnptions concerning pressure, viscosity, etc., and derive the coefficient 
values governing its evolution from the evolution observed. (This 'scout' process may be 
combined with thatfor diffusion coefficients: see below.) Scan data are not necessarily 
30 sufficient to imply values for aU dynamical variables. Motion being easier to observe at a 

distance than is pressure, it is usuaUy necessary here to specify boundary conditions imposed 
by the background and the injection device and solve an equation such as (1) or (3) together 
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with (2) above, for the tracked v and appropriate boundary conditions on pressure;?, as an 
equation in the less directly observed quantities/; and K. af (1) is replaced by an equation 
with more variables these are desirably included in this process.) For different injection 
plans, equation (1) can then be soh^ed to gnd v, as in the first approach to step (b). The 
visibility of the fluid to a scan may be a result of the detectabiHty of material dissolved or 
suspended in the fluid, if it is able to move fi-eely with the fluid (rather than be filtered, as 
would occur with large particles or with CO2 bubbles, the latter having an excellent visibility 
to ultrasound), or imposed by the scanning method itself such as the use in magnetic 
resonance imaging of preparing the hydrogen nucleus spins in some regions differently from 
those in others, so that the movement of the regions may be tracked. 

A third approach to step (b) exploits the real-time tracking of the material of final 
interest or the surrounding fluid, where naturally or by modification it is visible to the 
scamiing method used. This provides data from which the current v may be derived without 
assumptions concerning the determining physics. Unless the flow is constant, however, 
prediction of fixture values requires estimation of the equation coefBcients, as in the previous 
case. Such estimation involves solution of a boundaiy value problem fest enough to use the 
coefficient values in a bomidary-vahie computation of predicted v, which to be relevant must 
be in a time short compared to the injection process. This places higher demands on 
computing power. 

The prediction (c) rmy be performed by a variety of numerical solution methods, 
employing an already weU-developed art. . At the level of the continuum approximation' 
(where the finite number of diffiising particles is represented by an infinitely subdivisible 
range of concentration values, varying continuously from pomt to point) diflfiision of a 
molecular substance Xmay be governed by an equation 



^ = V . . VC) - V . (vC) + i [F, - /V, - R,^ (C, B)] 



(5) 



where D is flie diflRision tensor and the last term includes rate of addition of Xmolecules (if 
any), from blood vessels into brain tissue, rate of removal /v, of ^molecules by lymphatic 
dramage, and HnndfQB) is a kinetic term which accounts for the reaction of X molecules with 
compounds fixed in the tissue, such as binding receptors. These rates axe normaUy reported 
per unit volume of tissue; consequently we need the fector (1/^ in front of the last term in 
order to convert to a rate per unit fluid volume. If the substance JiT resists passage through tiie 
blood-brain barrier and does not react with the tissue, (5) reduces to 
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:^ = V.(2>«VC)-V.(vC) 

at ■ •" • 

or for difi&ision in static fluid, 

^=v.(D«ve) ^'^ 

St 

which we will use for illustration below. 

StUl another alternate type of model approximates the very large collection (order 10^^) of 
real difiusing particles by a much smaUer set of numericaUy simulated particles, each one 
moving probabiUsticaUy on a grid of points. In a time step ^, a particle's probabiUty of 
jumping fromapointtoaneighboring point ^ = x-.^,separatedfromit by the vector ^,.s 

proportional to exp(-/D-V2^0, y^^crc yHS -vAt) and D'' is the inverse of the matrixD. 
Figure 6 iUustrates this for the field of diffiision constants shown in Figure 5; a particle at 
point601 is more likely to jump, given a diffiisipn constant of8.0, than a particle at point 
602. Consequently, a set of particles injected at point 601 will spread out fester than a set at 
point602. The gradient operator VappUed to C in (ref.. static difiRise) corresponds to 
averagingthe net flow between^ and 3c (wherexcurrentlyhasmoreparticks than X, 
independentlyjumping,moreparticlesjumpfromxto x thanfi-om x tox),andthe 
divergence operator V. corresponds to summing the difference of jump rates into and out of 
X on opposite sides, in each direction. Loss and binding effects are included by terms that 
govern a particle's probability of disappearing or becoming bound or fi:eed;-infosion at given 
rate is modeled by the creation with controlled frequency of new particles to be tracked. 
Such particle models are a well-known field of study (see, for instance, Donald Greenspan, 
'Particle modeling', Birkhauser, Boston, 1997), and their solution may be summarized as 
"foUow the particles and see how many end up wha-e." 

The choice of a concentration-oriented model like (5) versus a probabilistic particle 
; model is dictated less by theoretical considerations than by available computing machinery. 
Any computer with pipelined vector architecture, such as the machines in the Cray senes. 
performs particularly fest on an array of concentrations, densities, etc., while a large system 
of less closely coupled processors operating in parallel may run festest by aUocating one or 

several particles to each. 
0 Either class of system can be optimized at the hardware level to create a computer 

dedicated to the fest solution of the particular class of problems arising in solution of 
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concentration problems and coefficient estimation of the lands necessarj' for this invention. 
The use of such hardware-optinized components in a clinical implementation fonns part of a 
preferred implementation, subject to economic feasibility. 

Coefficient estimation methods 

The first step in estimating coefficients is to exploit a priori knowledge of tissue 
types. For example, the high density of bone cortex, clearly apparent in CAT scans, is 
associated with eflfectively-zero diflSision and permeability tensors for the fluids and 
diffiising substances chiefly of interest here. These values may be assigned to the 
corresponding points in the grid of values for these tensors, with a high degree of confidence. 
This confidence level is itself recorded for these pomts, in an associated field of data that may 
be discarded once the estimate is complete. 

At other points, scan data are less determinative. For example, at the surfece of the 
cerebral cortex is a membrane with very low permeability, so that there is little diflusion 
across it. The membrane itselfis rarely visible in scan data such as MR images. In some 
parts of an image, its existence can be deduced by the presence of a boundary between 
cortical tissue and cerebrospinal fluid, so that a suitably tuned edge detection filter can with a 
high associated confidence level assign low permeability to surfece points. However, where 
two parts of the cortex are pressed together in a sulcus (fold), no local operator on ^ical 
current MR data can detect the existence of a boundary. The human observer may see that a 
membrane must be present, and approxhnately locate it, by continuation fi-om other parts of 
the image; but human input of such perception and reasoning is too labor-intensive for a 
clinical appKcation. As described in the Background, work exists on three-dimensional 
'vfrtual balloons', computed objects which fold themselves to fit the available data and 
stretch across such gaps m local evidence, with analogous results to human interpolation, but 
this is not yet routinely usable at clinical speeds. Where such techniques become available, 
implementations of the present invention may import data fi-om the resulting models of tissue 
stracture. but this is not part of the preferred first implementation. Other intelligence input 
may also be used, such as with touch screen fimctionahty, light wand fonctionality, or other 
30 user or automatic fimctionalities to add line features, surfece features, volume features 
composition (property) features to a model or even equations. 

In certain cases it will be essential to perform 'scout' injections, with difRising 
material whose concentration C(xuX2,xs.t) can be tracked over time in three dimensions. 



20 
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Taking as an example the static case. where.(7): can. be,used alone, its use requiring us to 
determinethe symmetric dififiision tensor , ;,vi .: . 



D = 



d^z 



22 ^^23 



/,3 dj2 djjj 



(8) 



in a region of interest (either spatially constant or varying from one grid point (y,k) to the 
next, as the phsyical situation requires). If Chas by notation the estimated value Qj^ at time /, 
then pointwise, the vector VC may be approximated by 

^ _r c -Crn), so that (7) gives us the discrete analogue 
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Catheters with multiple himens could be particularly use&l in this practice. A first scout 
injection could be released through one lumen and the tracking/difRjsion parameters 
determined. The actual active material may then be released through a second lumen (or the 
first lumen purged or emptied before administering the active material through the first 
lumen). 

For each time step we have vahies for Q*, and its neighbors, and hence for the vectors 
C^b , giving an instance where (7) is assumed approximately true. Given (for example only) 
SK such instances, if the C,^ are linearly independent, this suffices in principle to determine 
d and hence D uniquely. A similar process (requirmg more instances) can provide estimates 
of a larger set P of parameters, such as D and AT together. Where v is being estimated the 
algebra is largely similar, but can be separated fi-om the estimation of Z> because fluid motion 
leads to transfer of concentration which is antisymmetric between ne^bors (transfeiring 
from one point to another), where the difRision analysis above finds the symmetric influences 
on transfer. 

In practice, given the noisiness of data, it is better to find a least squares estimate (or 
other linear or exponential approximation) for D, using aU available instances of whose 
entries are non-negligible; and it is use&l to smooth the data C before entering this process. 
(The smoothing is best weighted by values fi-om an edge operator appKed to the background 
CAT. MR or other image data, since where there is evidence of sharp transition in the tissue, 
smoothing away sharp differences in Cis more likely to introduce error.) Other refinements' 
of this estimation process will be obvious to those skilled in the art. With digital imaging 
concurrent wife the procedure, pixel-based smoothing or blending are examples of modalities 
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for assessing or reflecting tissue transition. Edge feature sharpening techniques are more 
specific examples of such image enhancing treatment to reflect such tissue transitions. 

Where other solution methods for predicting C are preferred to direct manipulation of 
values C,jk, , or for other reasons as they arise, this estimation step may be performed by 
5 analogous but non-identical techniques, clear to those skilled in the art. For example, if a 
quantity g as described above is approximated by combinations 

c^Qi + C2Q2 + ^2^2 + • + ^i^^M of basis functions, estimation from data of coefficients in the 
equation used to predict the evolution of the a may give more robust results than estimation 
of pointwise values for parameters such as D. Such modifications fall within the scope of the 

10 invention here claimed, which primarily addresses the exploitation of any estimation method 
in creating a foil system for prediction and control. 

The condition above that the C,^ should be linearly independent is not a trivial one. 
Sfac randomly sampled vectors from a 6-dimensional space are likely to be independent, but 
the samples here are not drawn randomfy. In particular, if the medium is isotropic and the 

1 5 diflSising material spreads from a single point a, differences in simultaneous C value at two 
points xandy equidistant from a will typically represent only noise. Such data cannot 
directly reveal constants governing diflfiasion tangent to a sphere centered on a. However, the 
degeneracy of the estimation problem itself provides evidence of isotropy. For example, 
difiiision in an anisotropic medium generically has level sets of concentration that pass 

20 through a given point x in different directions at difSarent times. The estimation procedure 
can thus be reduced to fitting D as a scalar muifipfe of the identity matrix, for which the data 
will be sufBcient. 

Various forms of priori knowledge can strengthen this estimation process. For 
example, where (as with bone) a difi&ision tensor £> can be assigned with a high confidence 

25 value from the static image bacl^ound alone, a lower variance in the least squares fit to Z> 
should be required in modifying the estimate away from that assignment. In other case it 
may be known that a tissue type is anisotropic with characteristic diflfiisivities in high and low 
directions, but the direction of highest diflSisivity may not be deducible from the static 
backgroimd image. In this case, we would use the scout data (or incoming data as a plan is 

30 being executed) to estimate the orientation of £> at each point. Estimating only the three 
degrees of orientation freedom, with the eig^values assumed known, can be more robust. 
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Predictive control 

Image guidance permits a user injecting a substance to place a delivery instrument 
such as a catheter into (for example) a human brain, and see where it is located relative to 
scan-visible structures. In present practice, the user makes a plan for injection, determining 
5 the quantity of the substance to be loaded into the delivery system, the pressure or flow rate 
(perhaps variable) at wliidi the substance is to be iiyected, and the time over v/iuch the flow 
^ to be maintained. The objective of such a plan is to deliver the substance in desired 
quantities in or to 'target' tissues, often while minimizing the amount delivo^ to non-target 
tissues (where it would be wasted) and to vubi^ble tissues where the substance would do 

10 harm. It is also desirable to limit excess pressure, and the resulting edema (swelling) of 
tissue, in spati^ extent, intensity and duration. 

Commonly, the plan is ^tered into a computer tiiat will control the injection process, 
but direct hand control is also possible. This plan is then followed, with a change only if 
visually inspected images (e.g., U.S. Patent Nos. 6,026,316; and 5,964,705) of the diflRising 

15 substance makes clear that concentration, pressure or edema is not following the course 
expected by the user. Smce the evolving concentration is a 3D scalar field in the midst of 
complex 3D structures, visual inspection requires a strong grasp of the 3D relationships 
revealed by scan data. Currrait 3D display techniques do not display these relationships 
clearly enough to guarantee fest comprehension and appropriately swift action by the medical 

20 user. It is rare that undesired consequences are seen fest enough to limit their scope. 

The present invention replaces or augments such visual inspection by enabling the 
computer to enter this part of the. control loop. At the simplest level of control in the present 
invention, the user chooses a position or target area, and specifies an injection plan (quantity, 
duration, pressure/flow-rate) to the computer. In the novel process provided by the present 
.25 invention, the computer then solves or determines the transport properties on the assumption 
that this plan is followed, using the field of parameters P established as above together with 
available subsidiaiy data such as blood pressure to specify boundary conditions. The system 
displays the predicted values of concentration, pressure and edema, and the user has the 
opportunity to examine them at non-crisis speed and determme whether they are satisfectoiy. 

30 If not, the user repeatedly changes the device position (actually or virtually) and/or the 

specified plan, until an acceptable result is predicted. At this point, the plan is inplemented 
under computer control. 
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At a more strongly supported level of control, and the user specifies target tissues by 
marking them (for example, by using a 2D or mouse or other input deAdce such as a 
joystick to cUck on points, specifying ball-shaped regions of adjustable radius; or, outlining a 
target region; or, moving a 2D or 3D cursor tteough an image of the region of interest and 
5 seeing a resulting spray-can-like superposed image showing the intensity of target values as 
input for corresponding points; etc.). Minimum and maximum target concentrations, and 
limits on pressure, edema and their durations may be input separately. Similarly, the user 
may mark vulnerable tissues, and quantify their vulnerability. The user inputs a pian as 
above, and the computer predicts the resulting concentrations. In this mode, the conqjuter 
10 display includes markers for agreement with the user's goals, so that visual mspection is 
guided to possible regions of concern. When a chosen plan is implemented, using scan- 
visibility of the injected substance the computer can monitor deviation from the expected 
time course of concentration, and detecting problems earUer (and thus more usefally) than 
visual monitoring. Where deviation beyond a threshold level is detected, the computer 
15 derives revised estimates for the transport parameters P and deduces revised predictions for 
the results of the plan. If these involve values pre^specified by the user as unacceptable, the 
computer may pause the injection for new user input. If they are acceptable, the computer 
may then implement the procedure* ■ ; 

At a higher level again, where the computer is capable of carrying out multiple 
20 predictive simulations of transport in available time*' the computer may search the space of 
possible plans for the plan(s) with the result most quantitatively desirable according to the 
user-specified criteria. (For example, it may befgin with a heuristically chosen starting plan, 
perhaps a standard one, investigating the predicted result of successive changes in the plan, 
making those changes with improved results. In general the physics will not aUow exact 
25 achievement of a particular user-specified target concentration, so this is not well-posed as an 
inverse problem. The system's goal must be to find a plan whose results are within user- 
specified limits.) Such plan(s) may be presented to the user, to be accq)ted or rejected by 
key-strokes, mouse cKck, voice, or other input means devoid of quantitative detafl; or, the 
user may quantitatively modify the plan(s), view the consequences of each modified plan, 

30 ami select among this larger set. 

Where the computer is thus capable of generating plans adjusted to user goals, when 
monitoring the real injection reveals a deviation firom the expected concentration changes the 
computer may not only deduce corrections to the assumed parameter field P, create revised 
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predictions and provide a warning if the results are unacceptable, but find a revised plan with 
more satis&ctory results and (subject to user approval, or autonomously) implement the 
revisions while the injection is in progress. This variant of the invention requires the greatest 
speed, and is the most likely to require m^jlementation with dedicated hardware optimized 
S for the purpose of solving the necessary equations &st. 

Flow of the method 

In summary. Figure 7 shows a logic of the present invention. First, the target tissue is 
scanned (701), using one or more 3D imaging modalities such as CAT, ultrasound, fimctional 

10 magnetic resonance, etc. For more than one modality, the data may be brou^t (702) into a 
shared coordinate frame. Using the resulting information about tissue properties at each point 
in the scamied region, the points are classified into tissue types (703) and assigned transport 
parameter values according to the general properties of these types (704). In some cases this 
provides sufScient information about tran^ort properties in such tissue of the material to be 

15 injected. Commonly, the user will also perform scout injections and scan the resulting 
changes in concentration of the same or a similar material (710). This will thereby derive 
(71 1) estimates more specific to the patient, for the sul»-egions in which the scout injections 
provide robust data, as attested by low measures of noise. These two estimations of parameter 
fields are merged (715), each point's parameter values being more strongly influenced by the 

20 more confidently derived, estimate at, that point, and sharp changes being more acceptable 
where there is. evidence from any modality of a boundaiy. 

The planning specific to the, injection begins with choosing regions that should 
receive the injected substance, or that are vulnerable to it (718), and limits on allowable 
pressure and edema. This may be in the mind of the user, setting the mental stage for later 

25 decisions, or input to the computer by means of a graphical interfece. The user or in some 
implementations the computer then fixes a position for the instrument used (720), which at 
this stage may be either a physical image-guided process or a virtual placement within a 
displayed scan image. In the next step (721) the user or the computer specifies quantity, 
duration, pressure/flow-rate of the injection, and in (722) the computer predicts the resulting 

30 time course of concentration, using the transport parameters established at step (715). 

The user or the computer inspects the predicted values and determines (730) whether 
these are acceptable in the Kght of the injection goals. If 'No', the procedure returns to step 
(720), where the instrument position may be readjusted, and defines (721) new values for 
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quantity, duration, pressure/flow-rate of the injection. If 'Yes', the physical procedure begins: 
if the position was chosen virtuaHy, the instrament must first be placed there physicaUy 
(740). 

Step (740) initiates the injection proper, leading to (745) the first time unit (next after 

5 none) of the planned flow. As each time unit is completed a counter tests (750) whether the 
end of the plan has been reached. If it has, the injection procedure ends (780) with continued 
monitoring until the transport speed of the injected material has become negUgible. If not, the 
image guidance system monitors the concentration (751) and the computer compares the 
values with those predicted (760). If disagreement is within a pre-set tolerance, the procedure 

10 goes to the next step in the injection plan (745). If it exceeds this tolerance, the computer 
revises parameter estimates in the light of the new data (761), and makes new predictions 
with specifications of their confidence level. If these are unacceptable and there is not the 
time to prepare a satisfectory revised strategy, a test (762) diverts the procedure to end (780). 
If revision is practical, the user or the computer specifies a modification of the plan (771), 

15 computes the results (772) and tests (773) their acceptability. If they are satisfectory, the 
procedure continues (745) with the revised plan. If not, the system decides (774) whether 
time is running out. If it is, it goes to the emergency stop subprocess (762). If time is 
available, the search for an improved plan continues with a return to step (771). Thesejoops 
continue until the end process (780) is reached. 

20 The method ofdesignofthe release bfthe material fi-om the medical device may be 

based upon many different release structures or mechanisms. The release may be through 
actual holes or apertures in a cover for a hihieh or on the catheter, membranes on the.^heter, 
open/close stops or blocks on lumens in the catheter, and the like. Each of these different, 
types of systems will provide release or deUveiy characteristics that should be considered in 
25 the analysis ofthedeUvery rate profiles in the practice of the invention. 



BNSDCXJID- <WO ^0185230A2J_> 



wo 01/85230 
What is claimed is: 



21 



:'CT/US01/14633 



1 . A method for the interactive sunulation of the movement of a first material in an 
organism, oonprising: 

5 (a) constructing a field of local phj/sical quantities appropriate to motion, appearance 

and disappearance of said materkl in said organism, and 

(b) constructing a field of a rate at which fluid flowing m said organism as a 
consequence of bodUy processes or of an mjection process moves at eadi position, thereby 
transporting said material, and 
1 0 (c) constructing a model of a source fi-om ^vbach said material enters said organism, 

with any pressure differential or flux rate that causes it to enter said organism, and 

(d) computing a concentration field over time of said material, assuming a planned 
time course of said pressure differential or flux rate, and 

(e) using said concentration field to implement control of concentration of said * 
1 5 material by adjustment of the source fi-om which said material enters said organism. 

2. The method of Claim 1 wherein said movement that is simulated comprises at least 
one mechanism of movement selected fi-om the group consisting of diffiision, flow, chemical 
creation or destruction of molecules, and absorption. 

20 

3. A method according to Claims 1 or 2, where said organism is a patient or 
experimental subject, and said implemented control of the concentration optimizes resulting 
concentrations in different tissues according to target values therapeutically or scientifically 
chosen. 

4. The method according to Claim 3 wherein implemented control optimizes resulting 
concentrations at tissue targeted for therapy. 

5 . The method according to Claim 3 wherein implemented control comprises 
30 minimizing resulting concoitrations at tissue that is not targeted for therapy. 

6. The method according to Claim 4 wherein hnplemented control also comprises 
minimizing resulting concentrations at tissue that is not targeted for therapy. 
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7. A method according to Claim 1 where the field of local physical quantities is 
constructed from at least source selected from the group consisting ofc 

a) scan data that, by appropriate scaling or other physically caUbrated 
5 transformation, imply values for concentrations and/or motions of a 

material that has physical, chemical and biological properties different 
from those of the first material; 

b) magnetic resonance data that, by appropriate scaling or other physically 
caUbrated transformation, imply values for concentrations and/or motions 

10 ofa material that has physical, chemical and biological properties different 

from those of the first material to be predictive of its transport properties; 

c) by appropriate scaling or other physicaUy cahT>rated transformation, fix>m 
scan data recorded from a previous administration to the specific target 
region of that organism of material that has physical, chemical and 

15 biological properties different from those ofthe first material; 

d) scan data that permit segmentation of a region occupied by said organism 
into anatomical types, for which anatomical types of regions general - 

diffusion tensor values and absorption rates are known. 

20 8. A method according to Claim 2 where the field of difRision tensor values and of rates 
of absorption, creation and destruction is 6onstrdicted from at leat onesource of data selected 
fi-om the group consisting o£ 

a) scan data that, by appropriate scaling or other physiCaUy caUbrated 
transformation, imply values for concentrations and/or motions ofa 

25 material that has physical, chanical and biological properties different 

from those of the first material; 

b) rates of absorption, creation and destruction is constructed from magnetic 
resonance data that, by appropriate scaling or other physically caKbrated 
transformation, imply values for concentrations and/or motions ofa 

30 material that has physical, chemical and biological properties different 

from those of the first material; 

c) by appropriate scaling or other physically caUbrated transformation, from 
scan data recorded from a previous administration to the specific target 
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region of that organism of mataial that has physical, chemical and 
biological properties dififerent from those of the first material; or 
d) where the field of local physical quantities is constructed from scan data 
that permit segmentation of a region occupied by said organism mto 
anatomical types, for which anatomical types of regions general diflRision 
tensor values and absorption rates are known. 

9. The method according to Claim 7 wherein said general diflfiision tensor values and 
absorption rates are known from actual measurements on other organisms or wh&& the 
concentration field is computed by solving a finite differential equation constructed from a 
continuum model of some or aU of the processes of transport, absoiption, creation and 
destruction. 

9. The method according to Claim 7 wherein said general local physical quantities are 
known from actual measurements on other organisms or where the concentration field is 
computed by the use of a multiplicity of processors acting in paraUel, each handling a sub- 
region of the simulated volume, and using shared or separate memory of the field of local 
physical quantities. 



20 11. A method according to Claim 1 where the concentration field is computed by 
computing the motion of a number of simulated particles of said first material, with 
probabiKty of motion in a particular direction controfled by a fluid velocity field and a 
diflfiision tensor field and probability of absorption, creation or destniction controUed by 
other parts of the field of local physical quantities. 



25 
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12. A method according to Claim 1 1 where 

a) the concentration field is computed by the use of a multiplicity of 

processors acting in parallel, each processor handling a sub- 
population of one or more simulated particles, and at least some of 
said mult^licity of processors using shared or separate memory of 
the field of local physical quantities and where the multiplicity of 
processors is embedded on a dedicated chip or assembly of chips; or 
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b) where a computer explores a range of possible planned time courses 

for pressure diflfercntiaJ or flux rate and identifies certain options as 

more satisfectoiy than other planned time courses in terms of 
predetermined targets, and permits the user to implement an offered 
time course fi-om among said certain options by an act of assent or 
selection in a user interfece of software embodying the method. 



13. A method according to Claim 7 where comparison of a predicted 
concentration with real-time measurement of actual concentration is used to modify the field 
10 of local physical quantities to values and rates that predict a continuing change in 
concentration of said first material, and in pressure and edema. 

14. A method according to Claim 1 where target minimum or maximum concentrations 
and maximum pressures or extents of edema at particular locations are specified by the user, 
15 and a computer warns if said target minima or maxima are predicted to be unmet if a current 
planned time course of pressure differential or flux rate is executed or continued. 



15. A method according to Claim 7 where predicted time courses of concentrations under the 
20 assumptions of a current plan and of suggested options are graphicaUy displayed by the use 

of imagery. /. . ■ - ^ .s------ = 

16. The method according to Claim 15 wherein the user controls a region or slice that is 
displayed by the use of imagery or wherein the user controls optical properties relative to 

25 tissue regions and material concentrations are visualized. 

17. The method according to claims 1, 2, 7 or 8 wherein said material is a dissolved drug, 
dispersed drug, emulsified drug, a suspension of cells, or a suspension of nanodevices. 

30 1 8. A method for providing a plan for the introduction of a first material to the tissue of a 
patient comprising: 

providing an image of an interior vohmie of a patient, the image having more than one 
type of bodily composition within a field of view of the image. 
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assigning or determining values for at least two of said more than one type of bodify 
composition within said field of view, said values relating to at least one property 
selected from the group consisting of a) solvation ability, b) difSisivity, c) 
absorptivity, d) adsorptivity, and transmissivity for a class of confounds or a specific 
5 compound, 

selecting a position withm the field of view for introduction of a first material, 
selecting a rate of introduction of the first material at said position, 
determining a pattern of first material movement from said position based upon at 
least the rate of introduction and at least one of said properties, 
10 comparing said pattern of first material movement to objectives for providing said 

first materied to tissue within said patient, and 

modifying at least one of a) said position, b) said rate of introduction, c) length of 
time of first material introduction, and d) physical properties of said first material. 

15 
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